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With the goal in mind of deriving a method to compute quantum corrections for 
the real-time evolution in quantum field theory, we analyze the problem from the 
perspective of the Wigner function. We argue that this provides the most natural 
way to justify and extend the classical approximation. A simple proposal is pre- 
sented that can allow to give systematic quantum corrections to the evolution of 
expectation values and/or an estimate of the errors committed when using the clas- 
sical approximation. The method is applied to the case of a few degrees of freedom 
and compared with other methods and with the exact quantum results. An analy- 
sis of the dependence of the numerical effort involved as a function of the number 
of variables is given, which allow us to be optimistic about its applicability in a 
quantum field theoretical context. 
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I. INTRODUCTION 

Quantum Phenomena underlie most of Modern Physics. Alongside it brings in the neces- 
sity of dealing with complex functions and interference phenomena. In particular, determin- 
ing the time evolution of a quantum system is relevant for many areas of Physics. When a 
large number of degrees of freedom is involved, numerical methods based on the integration 
of the Schroedinger equation fail. This, however, is the typical situation in Quantum Field 
Theory (QFT). 

When computing expectation values of operators in the vacuum or in the equilibrium 
state, one can use the path-integral approach. Furthermore, in the computation of the 
Wick-rotated Green functions of the theory, the complex quantum weight of each trajectory 
is transformed into a positive-definite probability weight. This allows the use of efficient 
standard importance sampling techniques, such as the Metropolis algorithm or other Monte 
Carlo techniques. This information allows the extraction of the spectrum and other proper- 
ties of the theory. The same applies when studying Quantum Field Theory at equilibrium. 
However, even in this situation, there are important exceptions in which the weights are 
not positive definite and the standard Monte Carlo methods fail. This is often referred as 
the sign problem. One example occurs in certain Quantum Field theories, as QCD, at finite 
chemical potential. Many methods have been devised to obtain relevant information in this 
situation, meeting partial success [ij. However, it is generally accepted that, despite the 
efforts, no fully satisfactory solution has been found. This is particularly unwelcome, since 
full-proof predictions in certain areas of Physics, which are of great relevance and timeliness 
(such as Heavy Ion collisions) are lacking. 

When studying the quantum evolution away from equilibrium or from a initial state which 
is not an eigenstate of the Hamiltonian the situation is even more severe. In this context, 
path integral methods for time-dependent expectation values follow from the Schwinger- 
Keldysh formalism ^\^- This allows systematic perturbative calculations. However, non- 
perturbative phenomena are often crucial and we lack an efficient numerical computational 
method to deal with this situation (for a review see Ref. [^]). 

This problem arises in different areas of Physics such as Nuclear Physics, Quantum Chem- 
istry, Quantum Optics, etc. One of these areas is Cosmology, which actually triggered the 
interest of the present authors in the problem. Quantum fluctuations play a role at different 
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instances in the early Universe. One such case is at the inflationary era, by generating the 
density fluctuations which act as sources for the anisotropics of the cosmic background ra- 
diation and the formation of structures. Many authors argued that the fluctuations develop 

n n 

from quantum to classical, and can be treated as classical at late times [5[-pJ|. Another in- 
teresting epoch which depends crucially on the understanding of the quantum evolution of a 
quantum field theory, is that of preheating and reheating after inflation Properties 
of the present Universe, such as baryon number density, gravitational waves or magnetic field 
remnants, might depend upon this dynamics. Obtaining reliable estimates demand an ap- 
propriate treatment of the quantum field theory evolution from an initial state after inflation 
to the fully thermalized reheated Universe. To estimate these effects, several authors 15|- 
26| have employed the so-called classical approximation (For a somewhat different context 



see also Refs. 



271]- 28|). This consists on treating certain modes of the quantum fields as 



random fields with deterministic dynamics given by the classical equations of motion of the 
system. The randomness is imprinted in the initial conditions of the fields, often determined 
by a particularly simple initial quantum state. The authors who employ the classical ap- 
proximation, often present arguments to determine which modes can indeed be treated as 
classical and those that cannot, leading to different types of initial conditions in both cases. 
The last point is crucial since in Quantum Field Theory there are ultraviolet divergences. 
An appropriate treatment has to deal with them through Renormalization. 

The present work originated with the goal of determining the region of validity of the clas- 
sical approximation, estimating the size of the errors induced by it, and hopefully even ob- 
taining a way to go beyond it, by incorporating quantum corrections to the results. For that 
purpose we have started by analyzing simple quantum systems with very few degrees of free- 
dom whose quantum evolution can be obtained by numerical integration of the Schroedinger 
equation. The initial conditions and the type of Hamiltonians considered are inspired by 
the field theoretical and cosmological applications. Thus, we focus upon quartic interactions 
and gaussian or thermal initial conditions. 

In all cases we have compared the results obtained by the classical approximation with 
the exact quantum evolution. Furthermore, we have also studied the results obtained with 
different proposals for incorporating quantum corrections, limiting ourselves to those that 
have any hope of being applicable to the quantum field theoretical case. In particular, a 
very interesting technique is the two-particle irreducible effective action supplemented with 
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a certain truncation of the number of diagrams involved (2PI method) 29 1- 30|]. This 
truncation can be based on the loop expansion or on the 1/N expansion 3l|. The latter 
behaves in a more stable fashion and has been used in our analysis. We remark that the 
traditional Hartree method can be considered a particular case of the 2PI method 32|]- 34 1. 
In any case the method is limited to the determination of the quantum evolution of certain 
observables, such as the 2-point correlation function of the system (see Ref. jj] for a more 
complete list of early references on the subject). 

Another technique which has been recently applied in the context of Quantum Field 
Theory is the complex Langevin method SSj- 38| . This is based on complexifying the fields 
and studying the dynamics of the field trajectories induced by a Langevin equation in an 
additional Langevin-time variable with a purely imaginary drift term. Instabilities are of- 
ten found in the numerical integration of this equation, although authors have given sev- 
eral recipes to avoid them jssl. Furthermore, good results have been reported in certain 
cases j40|-j41[, so that we thought it was very interesting to apply it to our examples. Un- 
fortunately, we seem to be in a situation in which convergence to the right solution does not 
apply. This question deserves future study. 

In parallel to the tests explained before, we present a method to quantify and incorporate 
quantum effects based upon the Wigner function 42|]. This is a pseudo-distribution function, 
whose expectation values give us the matrix elements of Weyl-ordered products of operators 
in the quantum state of the system. The function is real, but not positive definite (see 



Ref. 



44| 



br an account of all its properties). The Wigner function satisfies an evolution 



equation 43|] in time, which determines the time- dependence of all these expectation values. 



One of the advantages of this method, is that it is particularly simple to see what is the 
meaning of the classical evolution and what is the nature of quantum corrections. This 
will be explained in the next section. This observation is not novel, and has led different 



researchers in different fields to focus on the Wigner 



unction and its evolution equation when 



451]- 481]. One example, is in Nuclear Physics 



attempting to describe quantum evolution 

where several authors |49| [s^ have developed methodologies which are very similar in spirit 
to our goal (see also [46[). Unfortunately, the detailed techniques seem hard to extend to a 
large number of degrees and thus to Quantum field theory. In our particular proposal we 
have dedicated some time to study the way in which the numerical effort involved scales 
with the number of degrees of freedom. A power-like growth is acceptable even if it involves 
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an enormous computational effort. Experience teaches us that the development of computer 
technology and algorithms will diminish the load in due time. An exponential growth is 
a killer. Our results presented below are promising and seem to allow the computation of 
quantum corrections in typical situations relevant for cosmological applications. This will 
be addressed in a future paper of the present authors. The present paper is to be considered 
a pilot study, in which we have the advantage of knowing what the exact quantum result is. 
The full field theoretical case demands a much higher computational cost, in addition with 
the necessity of dealing with issues as renormalization, as mentioned previously. 



For a more complete account of the literature we should mention the work of 51|, in 
which the authors advocate the use of the Wigner function evolution equation in Quantum 
Field Theory and give ideas on possible techniques to obtain an explicit calculation of the 
quantum corrections in this setting. 



We conclude this section by describing the lay-out of the paper. In the next section we 
review the definition of the Wigner function and derive its time-evolution equation. We also 
explain the meaning of the classical approximation in this context. The following section 
explains how, in trying to derive the quantum Liouville equation as a Fokker-Planck equation 
associated to a Langevin process, the non-positive definite character leads to pathological 
properties of the stochastic force. With this idea in mind, in Section HV] we present a method 
to relate the equation to a regular markovian process for which standard sample methods 
can be applied. This together with a certain coarse-graining approximation allows to set up 
a procedure to calculate systematic quantum corrections to the evolution in powers of h"^. In 
the following section, we study several simple cases for which we analyze the accuracy of the 
classical approximation, the truncated 2PI method and our proposal of the previous section. 
Special attention is paid to estimate the capacity of the method to deal with discretized 
lattice approximations to quantum field theory. In the concluding section, we summarize 
our results and discuss the advantages and limitations of our proposal. 
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II. THE WIGNER FUNCTION EVOLUTION EQUATION 

In Quantum Mechanics the expectation values of Weyl-ordered products of operators can 



be computed in terms of the Wigner function W{x,p,t) as follows |42[||44|: 

mw{Q,pm = I ^f{x,p)w{x,p,t) (1) 

where fw{Qi P) means a Weyl-ordered product of position and momentum operators, whose 
classical limit is f{x,p). For a pure state, given the wave function the expression of 

the Wigner function is 

W{x,p,t) = I ^^*{x + y/2,t)^{x-y/2,t)e'Py/'' (2) 

This can be extended to mixed states associated to a density matrix p{t) as follows: 

W{x,p,t) = I ^{x-y/2\p{t)\x + y/2)e^^y'^ (3) 

In both cases the Wigner function is real but not necessarily positive definite. 

The Wigner function satisfies an evolution equation in time which depends on the form 
of the potential. Here we will derive it for the case of Hamiltonian of the form 

^ = £ + V'(x) (4) 

The corresponding Schroedinger equation satisfied by the wave function is 

ido^{x, t) = - A^vp(3;, t) + ^*(^, t) (5) 

From this equation, using the relation Eq. [21 one can derive the equation followed by the 
Wigner function: 

+ i / f + f , t) ^{x - % t) e'py/\V{x + f ) - V{x - f )) (6) 
One can transform the right hand side in an obvious way and obtain 

/ 1 4 ^ ^ ('-'^ ^ l> - "'^ - !>)) - f ■ - f ■ 

Integrating by parts one can substitute ^ by —ip/h. On the other hand, the factors of y 
can be replaced by —ih^. We end up with the equation 

„ , X P dW(x,p,t) if, ih d , , ih d \ 
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This equation is well-known j43| and receives several names in the literature: Moyal equation 
or quantum Liouville equation. 

Let us now restrict to a potential of the form 

V{x) = ^x' + ^x' (8) 

The operator involving V can be expanded to give 

'8p 24 ^ ^8p 
Then, we are finally led to an equation of the form: 



abty(x.p.t) = -ii^^'''fP-'> +rM^"''<;-P-'> -^^'^'''"■P-'' (10) 

m 8x 8p 24 8p-^ 

Notice that the first two terms do not contain h. As a matter of fact, if we neglect the last 
term, the solution to this partial differential equation is very simple. It is given by 

Wo{xo{x,p,t),pQ{x,p,t)) 

where the functions Xq and po are obtained by running back in time to time zero the classical 
equations of motion from a point (x, p) in phase-space at time t. This is the so-called classical 
approximation to the quantum evolution. 

The last term contains all quantum effects and has dramatic consequences. In particular, 
the Wigner function is not guaranteed to remain positive at all times. Thus, computing 
expectation values with the Wigner function can be very unstable numerically, because it 
comes from a cancellation of both positive and negative terms which might be much larger 
than the overall sum. This is a typical sign problem, which might render difficult to compute 
quantum expectation values by probability methods. However, if we start at t = from a 
positive Wigner function it might take some time until the negative part contributes sizably, 
and expectation values can be determined with reasonable accuracy. 

The size of the last term is, in principle, small in macroscopic terms, being proportional 
to h?. However, this depends very much on the size of the third derivative of the Wigner 
function. This a time-dependent function, but it is clear that the initial distribution has 
an important effect on the accuracy of the classical approximation at initial times. This 
can be tested in our particular cases, since one can numerically integrate the Schroedinger 
equation. 
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If one focuses upon expectation values, the size of quantum effects and the errors com- 
mitted by numerical integration of the Moyal equation can be quite different. It is to be 
expected that the accuracy of expectation values is better for operators involving Q alone, 
than for those involving P. 

Let us restrict to a gaussian initial distribution, which appears naturally in many ap- 
plications. Assuming for simplicity factorization of the distribution in x and p, the initial 
Wigner function takes the form: 

For a pure state, the two standard deviations should be related as follows: 

h 

For a mixed state this condition is relaxed. For example, for the density matrix of a harmonic 
oscillator at thermal equilibrium, this product is given by 

h 



axCTp — 



2tanh(na;/3/2) 

This interpolates between h/2 at low temperatures and kT/u at high temperatures. 

Sticking to the pure state case, and given the scales involved in the problem, one can form 
dimensionless quantities which control the relative importance of quantum effects. The first 
one r is the usual one formed by taking the ratio of a classical quantity with the dimensions 
of action, divided by H. In our present case, this quantity is 

r = #^ (12) 

One expects smaller quantum effects for large values of r. However, there is another com- 
bination which seems more directly related to the size of the quantum term in the equation 
for the Wigner function. This is given by 

. = ^ (13) 

To get some insight into the structure of the Wigner function, we can study certain limits. 
For example, one can consider an ultra-quantum limit, in which we neglect the classical h- 
independent terms in the equation satisfied by the Wigner function. The equation can be 
integrated exactly in this case, and the result is a one-dimensional integral: 

W^^{x,p,t) = ,^ ,3,, e--'/(2-^) / dzeM^P^ - i>^Qz' - ^IvV^} (14) 
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Ultra quantum Wigner tunction 
saddle point approximation 




-10 -5 5 p 10 15 20 25 30 

FIG. 1. The Wigner function W^q in the ultraquantum hmit, for fixed values of x and t. 

where Q = —Xh'^xt/2A. This can be related to Airy functions. The shape of this function is 
displayed in Fig{T]for Q = 2 and 0"^ = 2. Notice the damped oscillations for positive p. It is 
clear that the Wigner function becomes negative in some regions, but the total integral is 
finite and positive. As a matter of fact, it is quite easy to understand this oscillatory pattern 
and its dependence on x and t, by evaluating the integral in the saddle point approximation. 
The result is also plotted in FiglH For small values of p, the approximation breaks down 
as expected, but it becomes quite precise for large values of which encompasses the 
oscillatory region. Introducing k, = 3Qp — crl/4:, the approximation is different for positive 
and negative values of k. In the first case we have 

2v^/€-i/4cos(2/€=^/V(27g2) - 7r/4) exp{~pal/{6Q) + ^^/(lOSg^)} (15) 

while for negative k, we have 

0F|«:|-i/^exp{-2|/s:|3/V(27Q2) -paJ/(6Q) + aj/(108g2)} (16) 

Notice that for large p the argument of the cosine is proportional to p^^'^/Q^^'^, so that it 
broadens for large times. Despite the complicated behavior of the Wigner function in this 
ultraquantum case, all expectation values of f{x) are time independent. 

One can go beyond this approximation by considering also the second term on the right- 
hand side of the Wigner function equation. This approximation is equivalent to the infinite 
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mass limit m — > oo. The new Wigner function is obtained from Wuq(x,p, t) by replacing p 
by p + V'{x)t. Again, all expectation values of x are time independent. 

MANY VARIABLES 

The previous results generalize to more than one variable in a fairly straightforward 
fashion. The equation satisfied by the Wigner function becomes: 



PndW dV dW\ y-^ d^V d^W 



° ' ' ^ \m dXn dXn dpn ) 24 ^ dXndXmdXr dpndpmdpr 



(17) 



Now we will consider two particularly important cases. The first one is a 0(N) symmetric 
potential: 

v-^\\x\\'+^mn' (18) 

With this potential the last term on the right-hand side of the Wigner function equation 
becomes 

n m ± 'I- ± m 

An interesting situation occurs when the initial distribution is 0(N) invariant. The Wigner 
function at all times would only depend on invariants {A = x ■ p, B = and C = 

The corresponding equation for W{x,p, t) = F{A, B, C, t) is given by 

- * (C^^S + ^AC^ + 4(2^^ + BC)^ + 8ABg + 

+ (2A^ + ^)c£k + (4^ + 8)^S) (20) 

If the initial Wigner function has typical values of A, B and C proportional to N, as in 
the case of independent variables, and we scale A to be proportional to 1/N, the quantum 
evolution preserves these properties. It is interesting to notice, that in this case the quantum- 
term in the evolution of the Wigner function is suppressed by one or two powers of N in 
the denominator. This means that in this particular large limit, the typical expansion 
parameters for the quantum evolution is 
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This is consistent with the conventional assertion that the large N dynamics is classical. 
Furthermore, notice that those terms containing third derivatives are suppressed by two 
powers of A^, instead of one. Thus, to leading order in the Wigner function satisfies a 
simplified equation containing only second derivatives. After some work one can write this 
leading quantum term as: 

The second case which we want to consider is that in which the coordinates are labeled 
by n, the points of a d-dimensional hypercubic lattice A. The coordinates will be referred 
as and the corresponding Hamiltonian is given by 



H 



The main property of this family of Hamiltonians is its invariance under the symmetry 
group of translations in d dimensions. Notice that this corresponds to the discretization 
of the Hamiltonian of a Xcf)^ scalar field theory in d-dimensions on a lattice of spacing a. 
The quantity 7r(n) is the conjugate momentum to 0(n), satisfying canonical commutation 
relations among them. If we apply the general formulas to derive the quantum Liouville 
equation for the Wigner function in this case, we obtain that the quantum term is given by: 

(23) 

Formally, it is possible to take the continuum limit a — y and write down the equation 
satisfied by the Wigner function in the case of quantum field theory. Here we will not be 
using it so we will not write it down explicitly. The reader can consult Ref. 5l| where it is 
spelled out. 



III. LANGEVIN APPROACH TO QUANTUM EVOLUTION 

In this section we will try to generate the Wigner function quantum evolution equation 
by means of a Langevin process. A typical Langevin dynamical process cannot do the job, 
since it preserves the positive character of the probability density. Hence, we need a non- 
trivial modification of the standard technique to apply to this case. In what follows we will 
see how this comes out exactly. 
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Let us begin by writing the classical evolution equations with the addition of a random 
force: 

i{t) - ^ (24) 
p{t)^-V'{x{t)) + F{t) (25) 

Our goal will be to study the properties of the force F[t) to recover the equation for the 
Wigner function. In order to do so in a simplified manner, let us discretize the time variable 
and write the equation relating x' = x{t + 5t) and p' = p{t + St) to x = x{t) and p = p{t): 

x' = x + Sx = x + dt^ (26) 
p' ^p + 5p^p-V\x')5t + 5F(t) (27) 

As wc will see, to spell out the x dependence of the random force, wc should write SF{t) = 
x^^^{t)r]{t), where the distribution of the variable rj{t) is controlled by a function p{r]). Now, 
let us write the distribution function for x' and p', which by definition is W{x',p';t + St). 
We get 

W{x',p';t + dt) ^ J drj J dx dp p{r])W{x,p,t) S{x' - x - dx) S{p' - p - Sp) (28) 

Now wc should eliminate the integrals over x and p with the use of the delta functions. For 
that purpose we should make use of the time-reversed evolution equations 

x^x'-5t^ = f{x',p',rj) (29) 

m 

p = p' + V'{x') 6t - x^'^ri = g{x\p\ r]) (30) 

We get 



W{x\p\t + 5t) = J dvpiv) Wifix',p',7]),g{x',p\v),t) J{x',p') (31) 

where J is the jacobian of the change of variables. Finally, we expand the equation in powers 
of St and T). Keeping terms linear in St only, we obtain a discretized version of the quantum 
Liouville equation provided we demand 

Jdr]p{rj) = 1 

{rf) = jdr}p{ri)rf^{) for < n 3 (32) 

and 

{ri')^ ldvp{v)v'-^ (33) 
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It is quite obvious that if the variable rj is real, the previous equations are incompatible with 
a positive definite distribution function p{r]), since in that case (77^") > or all moments 
should be 0. Despite this fact in the next section we will see how we can relate the equation to 
that of a truly positive-definite distribution function, for which sampling statistical methods 
are applicable. 

The generalization of the previous formulas to the case of several variables is fairly 
straightforward. In principle, the force is now replaced by a vector of forces SFi. These 
forces are functions of several random variables with non-positive definite distribution func- 
tions. A simple way to parametrize these forces is 

SFi = r)Xi (34) 

where the distribution function of p{r)) coincides with the one of a single variable. The re- 
maining variables Xi can be distributed according to standard positive-definite distributions. 

As an example consider the case of the 0(N) symmetric potential. In this case one can 
write Xi = ^JjT + and choose a simple positive definite distribution function p(r, ||.^||) to 
recover the time-discretized version of the Wigner-function equation. We leave the details 
to the reader. 

For the case of a d-dimensional lattice field potential given in the previous section, a choice 
like Xi — will do, provided the are independent random variables with vanishing 

average and {^f) = 1. 

IV. A NEW COMPUTATIONAL METHOD 

In the previous section we have seen how to reproduce the evolution equation for the 
Wigner function by means of a Langevin approach with a random force with non-positive 
definite distribution function. This is the starting point for a new procedure to approximate 
the quantum evolution which we will explain in this section. The method depends on three 
steps or approximations which are intimately connected among themselves. The first part 
is a coarse- grain approximation in the momenta pi, with a characteristic coarse-graining pa- 
rameter e. Next we will show how one can reproduce the effect of the non-positive definite 
distribution function p{r)) by means of a purely Markovian process involving ordinary prob- 
ability measures. This will allow the use of standard sampling techniques to generate the 
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distribution. The last step will be the introduction of a parameter k multiplying the quan- 
tum term in the evolution equation for the Wigner function. The parameter interpolates 
between a purely classical evolution (k = 0) and the full quantum evolution (for n = 1). 
The whole procedure can be used to compute the evolution of quantum expectation values 
in powers of k. This is similar to an expansion in powers of fi^, although part of the h 
dependence sits in the initial condition and is left unchanged. 

As mentioned in the previous paragraph, our first step is a coarse-grain approximation, 
which will amount to an approximation to the non-positive definite p{ri). For that purpose, 
we might relax the condition that (r^") = for n > 3. One possible realization of the 
conditions is achieved by the following family of distribution functions: 

Mm 

Pn{v) = S{v) + 2^ -^{^{V - eaO - ^iv + ea*)) (35) 

i=l 

where ji are positive numbers, aj are real values, and e is a free parameter. Although not 
explicitly indicated, the coefficients ■ji and do depend on N. They are determined by 
imposing that all even moments vanish and odd moments, given by 

Mm 

/^2p+i = {v"''-')n = 2 X: 7. ar'e'^^-'^ (36) 

i=l 

vanish for p < N, with the exception of p = 1 given by Eq. [331 If we take, without loss of 
generality, that the parameters are of order 1 , this condition implies that 

7i oc (37) 

In general, the solution to the set of constraints will not determine the parameters 7^ and 
ttj uniquely. This freedom in the choice of the parameters is a bonus, since one can test 
the effects of the coarse-graining on the results, by exploring different choices. Furthermore, 
a better approximation is obtained by taking larger values of A^, which will be referred as 
the level of the approximation. The number of terms M^r has to grow as A^ grows. An 
alternative method to improve the accuracy would be to reduce the value of e, thus reducing 
the effect of higher order derivatives of the Wigner function. As we will see later, there is 
a limitation to the minimal value of e, which is dictated by the range of time for which the 
method would be applicable. 

The next ingredient will be that of relating the Wigner function to a positive-definite 
distribution function, which can be approximated by samples. The same can be done for 
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the distribution function ^(77) as follows. Let us introduce a positive-definite function p{r], a) 
involving a discrete Ising-like variable a — ±1. This function is related to p{r]) by 



a=±l 

The function ^(77, a) is normalized as a probability distribution 

and hence the pref actor is given by: 

■^=<(T>p=^ J drip{ri,a)a 



(38) 



(39) 



(40) 



The construction can be extended to the coarse-grained version given before. Hence, we 
define 



The normalization condition implies 



l + a 



1=1 



(41) 



Mn 
i=l 



(42) 



A similar procedure can be applied to the Wigner function 



W{x, p; t) ^K{t)Y, W{x, p, a; t) a 



(43) 



(T=±l 



where W{x,p,a;t) is a well-defined probability distribution. This is equivalent to writing 
the Wi gner function as the difference of two positive definite functions. If we label by {)]^ 
the expectation values with respect to W{x,p, a; t), then the expectation values with respect 
to the Wigner function are given by 

{aO{x,p))^ 



dx dp 0{x,p) W{x,p; t) 



(44) 



w 



Now, one can obtain a time-discretized evolution equation for W{x,p, cr; t) involving p{r], a) 
as follows: 



W{x^p^a]t -\- 5t) = Y^ dr) p{r), p)W{x — 5x,p — 6p, p ■ a;t) 



(45) 
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where the displacements 5x-5p are those coming from the classical equations of motion with 
a force proportional r] (as before). Summing the previous equation over a and using the 
previous definitions, we re-obtain the discretized evolution equation for the Wigner function 
provided 

K{t + 5t) = MK{t) (46) 

We see that this implies that the normalization factor grows exponentially, and this is 
precisely the main numerical limitation of the method. 

Having set up an evolution equation for the probability distribution W{x,p, a; t) we can 
employ standard sampling techniques to approximate it. This leads us to the concept of 
signed samples: a collection of A4 points {x'"-* , } such that 

r 1 ^ 

dxdp J2 W{x,p,a;t)0{x,p;a) ^ —J20{x^'"\p^'''^;(r^''^) (47) 

a=±l a=l 

From Eq. HSjone can determine the time evolution of samples given a realization of the noise 
1]. This is a Markovian process where p(?7, /i) determines the conditional probability for a 
transition from one point in phase-space to a new one and a possible change of sign of the 
discrete Ising variable. 

If we use the coarse-grained distribution pjy^r],^), then with a certain probability we 
would simply evolve the system with the classical equations of motion, and with probability 
proportional to 6t we would produce a jump in the value of momentum and a possible flip 
of the sign of a. The reason for referring to the first approximation as coarse- graining in 
momentum space, has to do with the discrete magnitude of the jump (of order e) at each step 
of the time-evolution. If the Wigner function would be a polynomial in p, the approximation 
would be exact. 

If we start the evolution with a positive definite Wigner function, then ^-Y-W{x,p,t = 
0) = W{x,p,(7,t = 0). Hence, all points in the sample have cr*-"-' = 1. As time evolves some 
of the points acquire a negative value of a^""^ . At time t = nSt the probability that a point 
in the sample has negative a is given by 

P- = = ^ (1 - (< ^ >p)") - ^1 - ^'^'"''^'^^ (48) 

with A a number of order 1, which depends on the detailed form of pn- The probability 
approaches 1/2 exponentially in time. Once this happens we encounter a severe sign problem 
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in computing averages according to the formula Eq. HU since the averages result from strong 
cancellations from a = 1 and a = —1. At this point the method breaks down. The typical 
time when this happens is given by 

T.^ (49) 

Clearly things get worse as we decrease e and improve as we decrease A and h. However, 
the magnitude of e must be related to the typical range of variation in p of the Wigner 
function. Otherwise, corrections involving higher order derivatives of the Wigner function 
become large. The systematic errors associated to the coarse-graining can be checked by 
varying e or by using different choices of A^. For small enough e, the approximation should 
behave better for larger A^. 

Once the maximum acceptable value of e is selected, the method stated before only 
allows the computation of quantum expectation values for a range of times. This limitation, 
although important, does not destroy the usefulness of the method. Time-range limitations 
are already present in studying the classical evolution equation of a quantum field theory 
discretized on a lattice. The moment that the fluctuations become sizable at the cut-off scale, 
the discretized evolution equation fail to reproduce the continuum evolution. Fortunately, in 
many applications a good part of the interesting physical processes take place in a relatively 
short period of time. This is the case, for example, in the context of preheating in the early 
universe, which is one of the main motivations that we had for embarking in the present 
work. 

In situations in which a reliable full quantum evolution can only be carried for a too short 
lapse of time, the method can be used to compute first-order quantum corrections as follows. 
This is the third ingredient that we anticipated in the first paragraph of this section. The 
strategy is to consider a modified evolution equation with a new parameter k multiplying 
the last term of Eq. [TOl This can be interpreted as multiplying h by In this fashion one 
extends the time range of applicability of our method by the multiplicative factor l/n. Now, 
evolving the system for several other smaller values of k, one can determine the evolution 
of the expectation values as a power series in k. 

A practical problem concerns accuracy. Since, our goal is to estimate the quantum 
fluctuations, we realize that when reducing k, they have decreased by the same factor. To 
maintain the signal to noise ratio one should increase the size of the sample by /t^, which 
would increase the computation time by the same factor. Since the reduction in k was 
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dictated by the necessity to extend the range in time of the simulation, we conclude that 
this can be done at a cost in computer time which only grows polynomially with this range. 

Of course, the drawback is that one does not compute the full quantum effects but only 
to leading order in k (i.e. ^^). This by itself is an important result because it would serve as 
a measure of the size of quantum effects and as a next-to-leading correction to the classical 
approximation. Higher order powers of k are also computable at a higher cost in computer 
time. 

To test these ideas we decided to study several simple quantum mechanical systems 
for which the exact quantum evolution can be determined by numerical integration of the 
Schroedinger equation. The results will be presented in the next section and compared to 
other extensions of the classical approximation given in the literature. 

Before closing this section, we should comment on the modifications necessary to extend 
the procedure defined previously to the case of many degrees of freedom. Our method might 
not be optimal for the case of very few degrees of freedom, however, it was designed to be 
extensible in an affordable way to the case of many degrees of freedom. In this respect it 
differs from other proposals in the literature based on histogramming which seem impossible 
to extend to a large number of variables. 

The formal extension of all the approximations involved in the method to the case of 
several degrees of freedom is indeed trivial. In the way that the Langevin process was 
extended at the end of the last section, it turns out that there is a unique random variable 
rj having a non positive definite distribution function with similar or exact properties as the 
one appearing for the one- variable case. Thus, the coarse graining p — > and extension to 
positivity p{r]) — > p{r], a) remains exactly the same irrespective of the number of variables. 
The rest of random variables entering the force are of the conventional type and their number 
increases linearly with the number of degrees of freedom, and so does the computation time 
for a given time-step. The sample now, however, involves trajectories in the multidimensional 
state of the system. This is exactly the same as for the classical evolution (with random 
initial conditions), except that now there is a single additional Ising-like variable a. Thus, 
for a fixed number of trajectories the computational cost will grow in a similar fashion to 
that of the classical evolution. A new question to worry about is whether the number of 
trajectories needed to obtain results with a reasonable accuracy depend upon the number 
of degrees of freedom. This will be studied in the next chapter for one particular example. 
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V. TESTING THE CLASSICAL APPROXIMATION AND ITS EXTENSIONS 

In this section we will present the results of our tests of the classical approximation and 
of our proposed method to obtain quantum corrections, together with other proposals. The 
first cases are particularly simple situations with one or two degrees of freedom for which 
the exact quantum evolution can be obtained through the numerical integration of the 
Schroedinger equation. Later we will explore the first steps towards a possible application 
to quantum field theory. 

Our first example will be a simple anharmonic oscillator at intermediate values of the 
self-coupling. The potential is that given in Eq. [8] with the following choice of parameters: 

m = l ■ ^"^ = 0.5 ; A = 0.45 (50) 

The value of the dimensionless parameters mentioned in Section 2 are given by r = s = 1.57. 
We choose as initial condition a gaussian pure state with width given by (T/i2/3 = 0.45/21/3. 
Our main observable was taken to be the expectation value of the square of the position 
operator as a function of time, which is noted (Q^)(t). The Hamiltonian, the initial state 
and the observable were used in a previous paper with a similar spirit to ours 52I]. However, 
for illustrative purposes we are presenting the results for a higher value of the self-coupling 
A, for which quantum effects are stronger. 



The main results are collected in Fig. 2(a) 2(b) In the first figure the time evolution of 



the observable is displayed in units of the half-period of the A = system. This expectation 
value performs oscillations with a frequency close to that of the A = system. The classical 
approximation, also displayed in the figure, oscillates as well, but the amplitude gets damped 
very fast with time. This damping is a typical feature of the classical approximation which 
has been pointed out repeatedly (including Ref. 52]). The figure also shows two other curves. 
The first being the 2PI approximation obtained by keeping only the leading and next-to- 
leading (NLO) diagrams in a 1/N expansion jsil. Notice that in this case the amplitude of 
the oscillation is not decreasing, but there is a shift in the period oscillation. This might not 
be a serious drawback in extracting average properties over time. Finally, we also present 
the result of the new method explained before, which includes the classical approximation 
and the leading h"^ correction. The exact details are explained below. 



In Fig. 2(b) we present relative error of each approximation, namely the difference between 



the quantum evolution and the corresponding approximation, divided by the quantum result. 
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(a) 




(b) 



FIG. 2. (a): Time evolution of (Q^) compared to the classical approximation, the 2PI truncation 
to NLO order in 1/N expansion, and the method proposed in this paper, (b): Relative error 
committed in each of the approximations as a function of time. 



The first line corresponds to the classical approximation, which oscillates with increasing 
amplitude. Quantum corrections start being negligible and grow to a 20% level at 2t/T ^ 
0.8, and 40% level at 2t/T ^ 2.5. The 2PI curve appears to be the worst, but this is due 
to the shift in period with respect to the quantum curve. A more fair presentation should 
involve a comparison of the height of the maxima, for which 2PI is certainly better than 
the classical approximation. The last curve is our calculation including quantum effects 
up to linear order in h'^, using the method described in the previous section. Notice, that 
it provides a very accurate description of the quantum effects up to 2t/T ^ 1.4. Beyond 
this point it has more sizable errors, but certainly smaller than the classical approximation. 
Furthermore, it provides at least an estimate of the errors committed by employing the 
classical approximation. 

The actual procedure that we followed to determine the leading quantum correction 
(LQC) is the following. We re-scaled the size of the quantum term of quantum Liouville 
equation by using k. Then, we used the coarse-grained approximation to p{ri) up to the 
second level {N = 2, (rj^) = but (r/'') ^ 0), and the sampling method described in the 
previous chapter, to study the quantum evolution equation for a given value of k. The 
formula to obtain the approximation (LQC) including the contribution linear in h"^ to any 
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observable O is 

OlQC = Oelas + -(O, - Odas) (51) 

The curves depicted in Fig. |2]were obtained using k = 1/6. 

To check whether k = 1/6 is in the hnear regime, and to give an estimate of the size of 
the higher order terms in /i^, we repeated the procedure and obtained for several values 
of K (k = 0, 1/10, 1/8, 1/6, 1/5, 1/4, 1/2). In this way we get an idea about how the value 
of the observable interpolates between the classical (k = 0) and the quantum value = 1). 
In Fig |3] we display the result obtained for the expectation value of at the position of 
the third maximum (2t/T = 2.1). The y-axis gives the values obtained for the different 
values of n mentioned previously. We also display the value for the classical approximation 
(k = 0) and the full quantum result (k = 1). It is quite clear that the results follow an 
approximate linear dependence for k < 1/3. The straight line is the result of a linear fit 
(1 free parameter) in this range. The extrapolation to k = 1 of the straight line is very 
approximately our estimate of the quantum evolution up to next-to-leading order in (the 
leading order being the classical approximation). In this particular case we see that the 
linear quantum correction term (LQC) accounts for 80% of the quantum effects. Thus, with 
the addition of the classical approximation, we reproduce the exact quantum result with a 
3% error. In principle, one could go beyond the linear approximation and determine higher 
order corrections in fi^. If we add the result of k = 0.5 to the data and fit the results to a 
second degree polynomial in fi^ we get an even better approximation to the quantum result 
(second line in Fig. |3]). 

Although the previous results by themselves show that our procedure cannot be com- 
pletely misguided, we have analyzed the different sources of error in the determination of 
the correction presented above. Using jack-knife methods we can quantify the purely 
statistical errors. They increase with time but remain always at the level of a few percent. 
A much more difficult estimate is the effect of the coarse-graining in momenta. This can 
be estimated by changing the value of e and/or adding more terms in the discretization 
to impose {jf^'^^^) = 0. In particular, we have used results at e = 0.3 and three levels of 
discretization. The effect is more pronounced at the maxima and minima of the oscillations 
and the better the approximation, the closer the results to the actual quantum evolution. 
We estimate that, at most, errors (to the quantum correction) could be of the order of 
10-15% at the third maximum [2t/T = 2.1) rising up to 20-25% at the fourth maximum 
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FIG. 3. The value or {Q^) at the third maximum in time computed using several values of k, (see 
text). The lines are linear and quadratic fits. 



{2t/T = 3.0). The same conclusion follows both by comparison of the different levels as well 
as by extrapolation in (with a ~ 3) to e = 0. The conclusion is that, even if the errors 
are sizable, the method provides a good estimate of the quantum effects. 

Before embarking into the generalization to several variables, we tested the situation for 
another case having several distinct features which are present in some of the phenomeno- 
logical applications to cosmology. We considered a potential of the form 

V{x, t) = -^/i^ tanh(at)x^ + ^x^ (52) 

with the parameters chosen to be 

fj."^ = 2 ; a = 0.2 ; A = 0.4 (53) 

This time-dependent potential provides a smooth interpolation between a single-well and a 
double-well potential mimicking the situation occurring in hybrid inflationary models. No- 
tice that tunneling effects are possible now. The initial condition is a pure state given by a 



gaussian with o"yU^/^ = 1/4^/^. The same observable as before, (Q^), is displayed in Fig. 4(a) 



for a range of times for which the potential has basically evolved to the future asymptotic 
potential. Hence, as expected, the expectation value migrates from its initial value to per- 
forming oscillations around xl^, , where Xmin is the minimum of future asymptotic potential. 
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Also shown, is the corresponding curve for the same 2PI approximation mentioned earher. 
Although following the pattern of the quantum result, the differences are substantial. On 
the other, hand the classical evolution works quite well for this case. However, the addition 
of the Linear Quantum correction (LQC) using our method (with n = 1/2) makes the re- 



sult much better, as can be seen when looking at Fig. 4(b) , where we display the relative 
differences with respect to the quantum evolution as before. 

A full comparison of competitive methods suggested us to include the results of the 
Complex Langevin method described in 38|], and we invested considerable effort in doing 
so. The method generates a collection of histories as a function of an additional time variable. 
Our naive implementation, however, led to a growing number of trajectories blowing up in 
this additional time. It is easy to see that this is a feature of the complex evolution in the 
discretized new variable and in the absence of random noise. Both the problem and the 
possible cures have been documented in the literature of the subject. One can employ more 



refined discretizations or use a much smaller Langevin step [38l |. but this pays an obvious 
price in computational cost. A way out proposed in Ref. [39!] is to use an adaptive stepsize. 
Another possibility is a modification of the noise ^o|. Using these techniques we were able 
to eliminate most of the divergent trajectories. A small fraction, but growing in extra-time, 

38l | to discard them or go back in time 



remained. We took the attitude mentioned in Ref. 
and re-evolve them. 

Another problem mentioned in the literature is the issue of convergence. Sometimes the 
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system can show a limited decay to the equihbrium distribution or can converge to a wrong 



hmit 



40| 4l|]. This seems to be case in our tests of the previous examples. Furthermore, the 



results obtained seemed robust under changes of the methodology (adaptive step, different 
stepsizes, discarding or re-evolving) and stable under further evolution in the additional 
time. These results did not even show the right oscillatory pattern already seen in the 
classical approximation. Some authors 4l|] claim that to ensure the correct convergence 
other modifications are necessary. However, given that this was not the main issue of this 
paper, we decided to drop out these results and defer its study to further scrutiny. 



Extension to several variables 

Since our ultimate goal is that of studying the quantum evolution of fields, it is crucial 
to determine how does the new method that we have presented depend upon the number of 
degrees of freedom. The standard non-perturbative treatment of quantum fields proceeds 
through a lattice discretization and a subsequent continuum limit. Renormalization is a 
crucial ingredient in the process to obtain meaningful physical results. The latter aspect lies 
somewhat far from the scope of this paper and will be addressed in a future publication. The 
focus here is rather upon the numerical feasibility of the procedure. For attaining acceptable 
results one has to reach a number of variables within the range of those customary for these 
kind of simulations. A priori, the method presented here is capable of doing so, since its 
computational effort is comparable to that involved in the classical approximation for a 
given number of sampling trajectories, which has been used successfully in this context. 

However, there is a point of concern which we want to address. It might well happen 
that the number of trajectories needed to attain a given precision in the estimation of the 
quantum corrections grows with the number of degrees of freedom: A polynomial growth is 
acceptable, an exponential growth is not. 

As a testing example we have considered a lattice version of a two-dimensional scalar 



field theory which has been studied by other authors in this same context 53[. We take a 
real scalar quantum field in 1 + 1 dimensions with Hamiltonian 

H{t) = I dx[^n\x,t) + ^{mx,t)f + lfi'<P\x,t) + ^<P\x,t)] (54) 

where /^^ may depend on t. We have already explained the influence of h factors, so from 
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now on we assume natural units, h = c = 1. The dimensionless field and its conjugate 
momentum tt = (p satisfy equal-time canonical commutation relations 

[Tr{x), (f){y)] = -i 5 (x-y) (55) 

which gives Tr{x) dimensions of inverse length. 

The next step is to consider the lattice version of the previous Hamiltonian. Continuous 
space is approximated by a discrete number of points x„ = na separated by a distance a, 
the lattice spacing. To deal with a finite number of variables we must, in addition, put the 
system in a box of size L with periodic boundary conditions. Altogether, we end up having 
N = L/a variables (f)n{t) = (f){na,t). The corresponding conjugate momenta 7r„(t) satisfy 
the commutation relations 

[7rn,0m] = --^nm (56) 

where the factor of a is necessary to preserve the dimensions of the conjugate momentum. 
Naive discretization then leads to the Hamiltonian 

^"^ r1 1 1 A 



^ = « [2""" + 2 ^^'^"^ + 2-" + 24*^'] ^^^^ 

n=0 

with V0„ = {(pn+i ~ 4'n)/0'- After a suitable rescaling of the variables and of the parameters 
the Hamiltonian can be cast in the form Eq. [22l 

After presenting our system and its discretized version, let us consider the dynamical 
process that we will study following Ref. {ssl. The idea is to study the evolution of the 
system after a quench. In practice, this means that the /i^ parameter flips its sign abruptly 
at time t = from a positive value to a negative one. This can be seen as a limiting version 
of our previously smooth (tanh) transition from single to double-well. 

In practice, what we will consider is the evolution of the system for positive times starting 
(at t = 0) from an initial state given by the ground state of the Hamiltonian with A = and 
positive fi"^. This initial state is therefore gaussian as in previous examples, and is easily 
generated. For our numerical simulation we have taken the parameters of the model to be 

A = Sm^ a = 0.8/m (58) 



where the unit of mass m is given by m = Then we have studied this model for 

= 2,4,8,16 and 32. 
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For a real quantum field theory application one has to study the limit of a going to 0, 
with a suitable tuning of the parameters dictated by the renormalization conditions. For 
very small lattice spacings one might even encounter problems of critical slowing down, but 
these will affect other methods too. Here, we will be content with scaling the number of 
degrees of freedom and focusing on statistical significance and computational load alone. 

Finally, let us select our observables and present our results. Rather than working with 
the field variables we can work with their Fourier modes 0fc(t), since they decouple in the 
A = limit. We use the following normalization for the discrete Fourier transform 

N-l 

n' 

n=0 

where k = (27r/L)j for j = —N/2 + 1, —N/2 + 2, . . .., N/2. A similar expression applies for 
the modes Tik{t)- Reality of our original field implies = (p-kit) (and the same for tt). 

As mentioned previously for A = and positive fi"^ all the modes oscillate with a char- 
acteristic frequency ijj{k) = \/ Asiv? {jt^ / N) / a? + /i^. At the classical level, the flip in sign 
of /i^ produces that the low lying modes acquire an imaginary and start growing ex- 

ponentially. In the presence of a non-zero A the growth ceases once the non-linear effects 
become important. 

In Fig. |5] we present our results for the sum of expectation values of the square of each 
mode \(l)k\^{mt)) for two degrees of freedom N = 2. As a matter of fact this observable 
is just a discretized version of j dx (fP'{x). The time evolution as a function of mt is given, as 
obtained from the numerical integration of the N = 2 Schroedinger equation. The result of 
the classical approximation and of our LQC method obtained from k = 1/3 are also shown. 
The exact numerical integration has negligible errors at the scale of this and the following 
figures. The errors of the remaining approximations were obtained by applying a jack-knife 
method to the sample of trajectories. The total number of trajectories used for this data 
was = 8 X 10^ 

The results show the same pattern as before. The classical approximation captures the 
main features, but our LQC approximation calculation is capable of reducing the discrepancy 
substantially. Only at the latest times this difference exceeds the level of the statistical errors. 



In Fig. 6(a) and Fig. 6(b) we display the corresponding results for N = 8 and = 32 with 
a sample of size Ai = 4 x 10^. Here we do not have an exact result to compare with, so 
the main issue is the dependence of the errors on A^ for a fixed sample size. Errors of our 
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(b) 



FIG. 6. The same as Fig. [5] but for = 8 (a) and = 32 (b). 



method are larger than those of the classical approximation as expected, but do not seem 
to depend crucially on the number of variables. The intermediate values of A^, not shown, 
display exactly the same pattern. For all values of the relative difference between classical 
and linear QC approximation approaches a constant with errors diminishing as the sample 
size Ai grows. For the maximum values studied of order = 5 x 10'', we can estimate the 
quantum correction at our latest times ~ 3 with an accuracy of 10% without a significant 
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dependence on the number of degrees of freedom. 

In conclusion, the proposed method seems to scale reasonably well with the number 
of degrees of freedom. The computational cost is only a certain factor higher than that 
involved in the classical approximation. Thus, phenomenologically interesting applications 
are addressable within present high performance computing capabilities. 

VI. CONCLUSIONS AND OUTLOOK 

In this paper we have analyzed the real-time evolution of simple quantum systems. Both 
the form of the potential, as well as the type of initial conditions arc chosen to reflect relevant 
applications in Cosmology. The simplicity of the systems allows the numerical integration of 
the Schroedinger equation and, hence, can be used to test different approximation methods. 
The simplest one is given by the classical approximation which amounts to the classical 
evolution of a random variable distributed according to the initial wave-function. In our 
examples the classical approximation always performed well at initial times, capturing the 
qualitative features of the quantum evolution. When trying to improve on this approxima- 
tion it seems natural to focus on the Wigner function and the quantum Liouville evolution 
equation that it satisfies. It is simple to introduce a parameter k in the evolution equation 
that interpolates between the classical approximation {k — 0) and the full-quantum evolu- 
tion {k = 1). This parameter amounts to a rescaling in the value of fi^ appearing explicitly 
in the equation. 

We have presented a method based on samples and a discretization (coarse-graining) in 
the distribution in the conjugate momentum, and compared its results with the classical 
approximation, the next-to-leading 1/N truncation of the 2PI equations, and the numerical 
quantum evolution. The method can be used to determine the quantum corrections to 
the time-evolution of expectation values as a power series expansion in k. This provides a 
natural extension of the classical approximation (corresponding to the lowest- order). The 
results, for the examples considered, seem to capture a sizable part of the quantum effects, 
therefore providing a possible alternative to other approaches. Even when the method 
sizably departed from the exact quantum result, the discrepancy remained smaller and of 
the order of the quantum effect, and no anomalous instabilities were observed. Thus, it can 
at least serve to attach a level of precision to the classical approximation. 
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The method is certainly not a panacea. We emphasize that quantum averages are subject 
to the sign problem since the Wigner function is not positive definite. Our samphng method 
is equivalent to reweighting, which is certainly not a solution to the sign problem. However, 
if one starts with a positive definite distribution function, the problem only sets in at a 
later time, and reliable estimates of the quantum effects can be obtained at early stages of 
the evolution. The method can also be used to compute the full quantum effect by setting 
K — 1. However at fixed k, our coarse-grained sampling method severely breaks down 
beyond a critical range of times. However, one can go well beyond this point if one aims 
at computing the quantum corrections to order and not the full quantum evolution. In 
this case, the computational cost only grows in a power-like fashion with respect to the time 
range of the analysis. This can then be viewed as a systematic improvement with respect to 
the classical approximation, which at the least can serve to quantify the size of the quantum 
corrections involved. All other errors, arising from the statistical size of the sample or from 
the discretization in conjugate momenta are quantifiable. 

As emphasized in the introduction, our main goal is to be able to apply the method to 
the evolution of quantum fields in the early universe. For that purpose it is important to 
see how the computational cost depends on the number of degrees of freedom. Methods 
based on histogramming the Wigner function have costs that grow exponentially with the 
number of variables. In designing the methodology we focused on an approach which only 
grows in a power-like fashion, even if there are more efficient ways to handle systems with a 
few degrees of freedom. In the last section we have presented a pilot study to measure the 
rate at which the computational cost evolves with the number of variables. The comparison 
is done by monitoring the size of the sample in order to keep the degree of accuracy of the 
quantum corrections fixed. In our exploratory study we focused upon a two-dimensional 
quantum field theory case which has been subject of previous study as a toy model for 
cosmological apphcations. The model has been discretized to a system with degrees of 
freedom. We applied our technique to the model up to A^ = 32, and we found that the 
statistical errors on the quantum effects of our observables for a fixed sample size remain 
stable with the number of degrees of freedom. Of course, the simulation time grows with the 
number of degrees of freedom, as does for the case of the classical approximation. Even if 
computationally demanding, these calculations are feasible with present day technology, and 
so will be the case for our proposed method. The results presented in this paper allow us to 
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be optimistic enough to embark in a higher scale study with all the complexities involved in 
a quantum field theoretical treatment. 
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